function [CS_R_stack, CS_N_stack, gamma_hat_vec]=...
    Twostep_CS_calc(alpha,gamma_min_vec,data,theta_grid_array,Farray,giveMoments,giveWeight,giveVarianceEstim)
%Calculates two-step confidence sets. theta_grid_array
%is a cell array, where cell i contains the grid of values we want to
%consider to parameter i.  Here I assume we are interested in confidence
%sets for linear combinations of parameters, and the linear combinations of
%interest are given the Farray.  GiveMoments is used to pass the moment and
%derivate functions, while giveWeight is used to pass the desired weighting
%matrix and giveVarianceEstim passes a function used to calculate
%covariance matrices. alpha sets coverage, and gamma_min_vec sets the
%values gamma_min to be used in computing the two-step confidence sets
%(gamma_min_vec should be of the same dimension as Farray).

%Set number of simulations to use to critical values
crit_sims=10^6;


%Calculate dimension of moments
m=length(theta_grid_array);
theta_mat=theta_grid_array(1);
theta_mat=cell2mat(theta_mat)';
if m>1
    for n=2:m
        theta_old=theta_mat;
        theta_new=theta_grid_array(n);
        theta_new=cell2mat(theta_new)';
        theta_mat=[kron(ones(length(theta_new),1),theta_old) kron(theta_new,ones(length(theta_old),1))];
    end
end
moments=str2func(giveMoments);
gmat=moments(data,theta_mat(1,:)',0);
k=size(gmat,2);

%Calculate Wald confidence set
[CS_N_stack, Wald_store]=Wald_CS_calc(alpha,data,theta_grid_array,Farray,giveMoments,giveWeight,giveVarianceEstim);

%Calculate value of "a" (the linear combination weight) corresponding to
%gamma_min and different hypotheses
a_min_vec=zeros(length(Farray),1);
for n=1:length(Farray)
    F=Farray(:,:,n);
    F(:,sum(abs(F),1)==0)=[];
    p=size(F,2);
    
    s = RandStream('mt19937ar','Seed',1);
    RandStream.setGlobalStream(s)
    K_size=sum(randn(crit_sims,p).^2,2);
    J_size=sum(randn(crit_sims,k-p).^2,2);
    
    gamma_min=gamma_min_vec(n);
    
    crit_quant=@(a) quantile((1+a)*K_size+a*J_size,1-alpha-gamma_min);
    crit_obj=@(a) abs(chi2inv(1-alpha,p)-crit_quant(a));
    a_min_vec(n)=fminsearch(crit_obj,0);
end

[CS_R_stack, K_store, S_store, theta_mat]=Robust_CS_calc(alpha,a_min_vec,data,theta_grid_array,Farray,giveMoments,giveWeight,giveVarianceEstim);

%Calculate gamma_hat
gamma_hat_vec=zeros(length(Farray),1);
for n=1:length(Farray)
    F=Farray(:,:,n);
    F(:,sum(abs(F),1)==0)=[];
    p=size(F,2);
    
    s = RandStream('mt19937ar','Seed',1);
    RandStream.setGlobalStream(s)
    K_size=sum(randn(crit_sims,p).^2,2);
    J_size=sum(randn(crit_sims,k-p).^2,2);
    
    a_diff=(chi2inv(1-alpha,p)-K_store(:,n))./S_store.*(Wald_store(:,n)>chi2inv(1-alpha,p));
    a=max(a_diff);
    a_vec(n)=a;
    temp=mean((1+a)*K_size+a*J_size<chi2inv(1-alpha,p));
    gamma_tilde=1-alpha-temp;
    gamma_hat_vec(n)=max(gamma_tilde,gamma_min_vec(n));
end

end

